PLh 



Numerical models of rotating accretion flows around black 
holes 

Igor V. Igumenshchev^ 

-- -^ Astronomy & Astrophysics, Goteborg University and Chalmers 

^ . University of Technology, 41^ 96 Goteborg, Sweden 

0^ 

Abstract. Numerical, two-dimensional, time-dependent hydrodynami- 
^^ . cal models of geometrically thick accretion discs around black holes are 

^ ' presented. Accretion flows with non-effective radiation cooling (ADAFs) 

can be both convectively stable or unstable depending on the value of the 
OO ' viscosity parameter a. The high viscosity flows {a ~ 1) are stable and 

have a strong equatorial inflow and bipolar outflows. The low viscosity 
flows {a ^ 0.1) are convectively unstable and this induces quasi-periodic 
^ . variability. 
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Q> ■ Advection-dominated accretion flows (ADAFs) have recently attracted much 

(— I I attention since they naturally explain the properties of X-ray transients, low 

Q^' luminosity active galactic nuclei and other high energy objects. Most of the 

information of ADAF structure derives from a simplified vertically integrated 
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1. Introduction 



^ • approach, which reduces the complicated three-dimensional problem of accre- 

j/3 I tion flow hydrodynamics to a one-dimensional problem (see Narayan & Yi 1995; 

C^ ■ Chen, Abramowicz & Lasota 1997; Narayan, Kato & Honma 1997; Igumen- 

shchev, Abramowicz & Novikov 1998). In the vertically integrated approach only 
the radial structure of the disc is studied in a detailed way. Due to significant 
simplifications introduced by the vertical integration, some important effects, 

C^ ' such as convection (Narayan &; Yi 1994; Igumenshchev, Chen & Abramowicz 

1996) or outflows (see discussion in Narayan & Yi 1994) are not properly treated. 
Understanding of ADAF properties could be improved by a discussion of two- 
dimensional time-dependent hydro dynamical models, where one explicitly treats 
both the radial and vertical structure of the accretion flow. 

In this contribution we present preliminary results of two-dimensional nu- 
merical simulations of rotating accretion flows around black holes. Complete 
discussion of the results will be published (Igumenshchev & Abramowicz 1999). 
We construct hydrodynamical time-dependent models with viscosity pa- 
rameter a ~ 0.1 — 1 and high geometrical thickness. We assume a simplified 
model of the radiative cooling: (1 — e) of the energy generated by the viscous 
dissipation is radiated away. We consider a large value of e (0.5 < e < 1), when 
the accretion flow has a low efficiency of the conversion of its internal energy 
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to the escaping radiation. We demonstrate that stability of the flow strongly 
depends on the value of viscosity. The low viscosity flows (a ^ 0.1) are convec- 
tively unstable, and the instability produces a quasi-periodic behaviour of the 
accretion flows and outflows. In the case of high viscosity (a ~ 1), the convective 
instability is suppressed, and the flow is stable. 

2. Numerical method 

Our numerical technique is based on the solution of the non-relativistic Navier- 
Stokes equations in the spherical coordinate system {r, 9, cp) in a stationary 
gravitational field of the black hole. The gravity of the black hole is modeled by 
the Newtonian potential 

^W = -|^, (1) 

2 r 

where Vg = 2GM/c^ is the gravitational radius of black hole of mass M. We 
assume axial symmetry of the flow with respect to the rotational axis that co- 
incides with the = direction. We take into account the contribution of 
all components of the viscous stress tensor to the equations of motion and the 
energy equation. Shear viscosity is only considered. 

For simplicity, we describe the dynamics of accreting plasma in the frame- 
work of the one fluid approximation. We use the equation of state for an ideal 
gas. The adiabatic index of gas is assumed to be a constant, and takes the value 
7 = 3/2. In this case the equation of thermal energy conservation can be written 
in the form: 

7 O 

where the operation d/dt is the comoving (Lagrangian) time derivative, S is the 
specific entropy, Qmsc is the dissipation function and Qrad is the volume cooling 
rate. The problem of the radiation losses in the high temperature magnetized 
plasma is a quite difficult one and it is far from the complete solution. We do not 
address this problem here, but instead assume a simple model for the cooling 
rate, 

where e is a parameter, e < 1. The case e = 1 corresponds to the non-radiating 
accretion flow, whereas e = corresponds to the isentropic flow. 
The kinematic viscosity coefficient is described by 

u = a-^, (4) 



where a is a constant, < a ^ 1, c^ = \JRT / ^i is the isothermal sound speed. 



and Q.K = cJrgl^r"^ is the Keplerian angular velocity. 

To solve the Navier-Stokes equations we split the numerical procedure of 
one time step At calculation into two parts, hydrodynamical and viscous. The 
hydrodynamical part is calculated by using the explicit finite-difference PPM 
algorithm developed by Colella & Woodward (1984). The viscous part is solved 



by the implicit operator splitting method for the contributions to the equations 
of motion. The viscous contribution to the thermal energy equation (2) is cal- 
culated using an explicit scheme. The time step Ai is chosen in accordance with 
the Courant condition for the hydrodynamical sub-step. 

We use the absorbing inner boundary at the radius 3rg , which is the location 
of the last stable orbit of the Schwarzschild black hole. The outer boundary is 
located at the radius ~ SOOr^. At this radius we inject the matter with a fixed 
rate Mj„j, close to the equatorial plane. Injected matter has zero r- and 9- 
components of velocities, and the angular momentum close to the Keplerian 
one. 

We look for stationary or quasi-stationary solutions, which are obtained in 
the course of the time-dependent calculations from the initial state. This requires 
to follow the evolution of flow pattern during a few characteristic accretion times, 
measured as an average time of motion of the matter from the outer boundary 
to the inner one. 

3. Results 

We have calculated a number of evolved models with four different values of 
a = 1, 0.3, 0.1 and 0.03, and different values of e, which vary in range 0.5 — 1. 
The models show a weak dependence on e, and they strongly depend on a. 

For large value of a = 1 and 0.3 the models are stable. They are symmetric 
with respect to the equatorial plane and show no time-dependent behaviour. 
Figure || represents the flow pattern of the model with a = 1 and e = 1. In 
this figure the meridional cross-section of the model is shown, and vertical axis 
coincides with the axis of rotation. Black hole locates in the origin (0,0). Axes 
are labeled in the units of rg. Upper left plot shows the contours of density p. 
The contour lines are spaced with Alogp = 0.1. Upper right plot shows the 
contours of pressure P. The lines are spaced with AlogP = 0.1. The density 
and pressure monotonically increase toward the black hole. In lower left plot 
the arrows with the length in relative units show the momentum vectors pv. 
The flow pattern consists of the equatorial inflow and bipolar outflows, which 
originate very close to the black hole, at radius 8rg. Lower right plot shows 
the contours of Mach number Ai. The lines are spaced with AA4 = 0.05. The 
maximum value of A4 at given radius is reached at the equatorial plane. Two 
stagnation points, where A^ = 0, locate at the axis of rotation at radius 8rg. 
The flow is subsonic {A4 < 1) everywhere in our computation domain, which 
has inner boundary at r = 3rg. 

The model with a = 0.3 does not show deep outflows. The stagnation points 
locate at radius 80 — 90rg on the axis of rotation. Inside this radius the matter 
moves to the black hole almost radially. Oppositely to the case of a = 1 the 
distribution of the Mach number in this model has an equatorial minimum at a 
given radius. This minimum indicates that the inflowing matter is overheated in 
the inner equatorial parts due to viscous dissipation. In the low viscosity models 
this overheating is the reason for development of the convective instability. 

Low viscosity models with a < 0.1 are unstable. They demonstrate rich 
and complicated time-dependent variations of flow pattern. Example snapshot 
of flow pattern of the convectively unstable model with a = 0.1 and e = 1 is 
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Figure 1. The flow pattern of stable model with a = 1 and e 
(see text for details). 




I t 



\\\\ 1 1 1 1 1 ' ' 




Figure 2. The flow pattern of convectively unstable model with a 
0.1 and E = \ (see text for details). 



shown in Figure ^. In the upper left plot for density distribution it is clearly seen 
the reason for the instability — hot convective bubbles, which quasi-periodically 
originate in the innermost region (4 — Gr^ from the center) of accretion flow and 
propagate outward. 

The bubbles arise from the initial perturbations, which usually appear in 
the ^-directions, where the maximum mass flux onto the black hole occurs at 
that moment. As a rule, the bubbles originate near, slightly above or below, 
the equatorial plane. When a bubble has developed, it forces the direction of 
the maximum mass inflow to change. In the new direction, a new convective 
bubble originates when the previous bubble reaches the radial distance 50-100 Vg. 
This cycle repeats quasi-periodically. Typically, convective instability produces 
sequences of convective bubbles, where the previous and next bubbles originate 
and move outward in different hemispheres with respect to the equatorial plane. 
The quasi-periodic behaviour of the convective bubbles, accompanied by the 



outflows in polar directions, produces a signiflcant variability of the flow pattern 
inside r ~ lOOr^. The intensity and directions (up or down) of the outflows 
strongly correlate with the convective activity in the vicinity of the black hole. 

The spatial scale of perturbed motion in the accretion flow becomes smaller 
with decreasing of the viscosity parameter a. Accordingly, the flow pattern is 
more complicated in model with a = 0.03. Small scale vortices accompany the 
convective motion. The vortices exist quite a long time and strongly interacts 
with convective bubbles. 

We have studied the time-dependent behaviour of the unstable models cal- 
culating radiative energy losses in the case of proton-electron bremsstrahlung 
cooling for optically thin plasma, 

L{t) (X f p^T^/^dV. (5) 

Here the integration is taken over the volume of a sphere with the radius 100 r^. 
Fourier analysis of time series of the radiative cooling rates shows the presence 
of a strong feature at frequencies 10 — 20 (Mq/M) Hz in the power spectrum 
of model with a = 0.1. In the case of a = 0.03, a weaker feature is observed 
at 1 — 2 (Mq/M) Hz. These oscillations can be explained by quasi-periodic 
presence of convective bubbles and perturbations of the accretion flow introduced 
by these bubbles. Note, this variability has typical time scales in the observed 
quasi-periodic oscillations (QPOs) range of the Galactic black hole candidate 
X-ray sources. 
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